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Method for Compensating Mild Lateral Velocity Variations 
in Pre-Stack Time Migration in the Frequencv-Wavenumber Domain 

[0001] This application claims the benefit of U.S. Provisional Application 
No. 60/415,382 filed on October 2, 2002. 

FIELD OF THE INVENTION 

[0002] This invention relates generally to the field of seismic prospecting 
and, more particularly, to imaging of seismic data. Specifically, the invention 
is a method for prestack time migration of seismic data where the velocity 
function has mild lateral variation. 

BACKGROUND 

[0003] In the oil and gas industry, seismic prospecting techniques are 
commonly used to aid in the search for and evaluation of subterranean 
hydrocarbon deposits. A seismic prospecting operation consists of three 
separate stages: data acquisition, data processing, and data interpretation. 
The success of a seismic prospecting operation is dependent on satisfactory 
completion of all three stages. Petroleum engineers know how to produce 
hydrocarbons from reserves found by a successful prospecting operation. 

[0004] In the data acquisition stage, a seismic source is commonly used to 
generate a physical impulse known as a "seismic signal" that propagates into 
the earth and is at least partially reflected by subsurface seismic reflectors 
(i.e., interfaces between underground formations having different acoustic 
impedances). The reflected signals (known as "seismic reflections") are 
detected and recorded by an array of seismic receivers located at or near the 
surface of the earth, in an overlying body of water, or at known depths in 
boreholes. The seismic energy recorded by each seismic receiver is known 
as a "seismic data trace." 
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[0005] During the data processing stage, the raw seismic data traces 
recorded in the data acquisition stage are refined and enhanced using a 
variety of procedures that depend on the nature of the geologic structure 
being investigated and on the characteristics of the raw data traces 
themselves. In general, the purpose of the data processing stage is to 
produce an image of the subsurface geologic structure from the recorded 
seismic data for use during the data interpretation stage. The image is 
developed using theoretical and empirical models of the manner in which the 
seismic signals are transmitted into the earth, attenuated by the subsurface 
strata, and reflected from the geologic structures. The quality of the final 
product of the data processing stage is heavily dependent on the accuracy of 
the procedures used to process the data. 

[0006] The purpose of the data interpretation stage is to determine 
information about the subsurface geology of the earth from the processed 
seismic data. For example, data interpretation may be used to determine the 
general geologic structure of a subsurface region, or to locate potential 
hydrocarbon reservoirs, or to guide the development of an already discovered 
reservoir. Obviously, the data interpretation stage cannot be successful 
unless the processed seismic data provide an accurate representation of the 
subsurface geology. 

[0007] Typically, some form of seismic migration (also known as "imaging") 
must be performed during the data processing stage in order to accurately 
position the subsurface seismic reflectors. The need for seismic migration 
arises because variable seismic velocities and dipping reflectors cause 
seismic reflections in unmigrated seismic images to appear at incorrect 
locations. Seismic migration is an inversion operation in which the seismic 
reflections are moved or "migrated" to their true subsurface positions. 

[0008] There are many different seismic migration techniques. Some of 
these migration techniques are applied after common-midpoint (CMP) 
stacking of the data traces. (As is well known, CMP stacking is a data 
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processing procedure in which a plurality of seismic data traces having the 
same source-receiver midpoint but different offsets are summed to form a 
stacked data trace that simulates a zero-offset data trace for the midpoint in 
question.) Such "poststack" migration can be done, for example, by 
integration along diffraction curves (known as "KirchhofT' migration), by 
numerical finite difference or phase-shift downward-continuation of the 
wavefield, or by equivalent operations in frequency-wavenumber or other 
domains. 

[0009] Conversely, other seismic migration techniques are applied before 
stacking of the seismic data traces. In other words, these "prestack" migration 
techniques are applied to the individual nonzero-offset data traces and the 
migrated results are then stacked to form the final image. Prestack migration 
typically produces better images than poststack migration. However, prestack 
migration is generally much more expensive than poststack migration. 
Accordingly, the use of prestack migration has typically been limited to 
situations where poststack migration does not provide an acceptable result, 
e.g., where the reflectors are steeply dipping. 

[0010] In some cases, reflector dip can exceed 90 degrees. As is well 
known in the seismic prospecting art, it may be possible to image these 
"overturned" reflectors using data from seismic "turning rays." Prestack 
migration techniques must be used in order to obtain an accurate image of 
overturned reflectors from seismic turning ray data. 

[0011] There are two general types of prestack migration, prestack time 
migration and prestack depth migration. A background seismic wave 
propagation velocity model to describe the seismic wave propagation velocity 
in the subsurface is needed in the seismic imaging. In a region where the 
subsurface seismic wave velocity varies only in the vertical direction, the 
seismic imaging method used is pre-stack time migration (PSTM). In a region 
where the subsurface seismic wave propagation velocity varies in both 
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vertical and lateral (or horizontal) direction, pre-stack depth migration (PSDM), 
needs to be used to give accurate results. 

[0012] Both vertical time and lateral position errors occur in time migration 
whenever there are lateral velocity variations in the subsurface. See, for 
example, "Plumes: Response of Time Migration to Lateral Velocity Variation", 
by Bevc, et al.. Geophysics 60, pp. 1118-1127 (1995); or "Systematics of 
Time Migration Errors", by Black and Brzostowski, Geophysics 59, pp. 1419- 
1434 (1994). Mild lateral velocity variations may be approximately 
compensated by interpolating multiple velocity functions or multiple travel time 
tables from different locations when the PSTM is performed in space domain, 
as in the Kirchhoff summation algorithm. See, for example. Seismic Data 
Processing, by O. Yilmaz, Society of Exploration Geophysics, pages 269-271 
(1987). A common offset seismic section is defined as the seismic section of 
the same source-receiver spacing (offset). PSTM of a common offset section 
in the frequency-wavenumber(A:,ty) domain is much more efficient than in real 
space. See, for example, "Recursive Wavenumber-Frequency Migration", by 
Kim, et al.. Geophysics 54, pp. 319-329 (1989); or the Finn and Winbow 
reference cited below. However, there is no known way to treat mild lateral 
velocity variations in the Fourier domain similar to the methods used in the 
space domain. Very often, the water bottom may have a small dip in a marine 
survey, or the subsurface velocity may have an increasing or decreasing trend 
along a specific lateral direction. These lateral velocity variations are usually 
big enough to cause errors, particularly in the lateral position of steeply 
dipping reflectors, when the traditional PSTM is applied using a single 
laterally-invariant velocity function. Full PSDM will properly position seismic 
images. However it is far more expensive than PSTM. 

[0013] A current method used to compensate the errors caused by the 
mild lateral velocity variations is to do residual move-out correction (lining up 
images of different offsets by shifting the image vertically) before the final 
stack (adding images of different offsets together). The purpose of this 
residual move-out correction is to enhance the final stack image. It can not 



Express Mail Label No. EU959479572US 



Filed on July 31, 2003 



Attorney Docket No. PM 99.017 



correct the lateral position error. It can not focus diffraction energy correctly. 
Therefore, it will loose the sharpness of fault images and give incorrect 
reflector dips (Bevc, 1995). The correction to the lateral velocity variations 
should be incorporated into the migration operator in order to correct both the 
vertical time error and the lateral positioning error. 

[0014] What is needed is a practical method for accounting for mild lateral 
velocity variations in the frequency-wavenumber domain. The present 
invention provides such a method. 

SUMMARY OF THE INVENTION 

[0015] In one embodiment, the present invention is a method for pre-stack 
migration of common offset seismic traces obtained from a subterranean 
region where the seismic velocity varies both laterally (with x) and vertically 
(with said method comprising steps of (a) selecting seismic velocity 
functions v(^) at at least two lateral (jc^) locations (for a two-dimensional case 
and three lateral locations for a three-dimensional case) in said subterranean 
region, the number of such locations being determined by the degree of 
lateral velocity variation; (b) transforming said common offset seismic data 
traces from the space-time (x-t) domain to the wave number-frequency (k-co) 
domain; (c) calculating a travel time map for each Xs location using the velocity 
function selected for such location; (d) calculating for each travel time map 
from step (c) a corresponding map of r as a function of wave number p, 
where cot is the phase shift in the ^-(y domain corresponding to the migration 
time shift in the x-t domain; (e) using the r -maps to find t (p) as a linear 
function of x with a certain slope and intercept at each depth (^) location in 
the subterranean region and each value of p] (f) forming the migrated image 
from said seismic traces in the k-w domain using the r intercept and slope 
values from step (e) in a pre-stack time migration equation with k shifted by an 
amount equal to the r slope from step (e) multiplied by ay, and (g) reverse 
transforming said migrated image back to the space-time (x-t) domain. 
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DESCRIPTION OF THE DRAWINGS 

[0016] The patent or application file cx)ntains at least one drawing 
executed in color. Copies of this patent or patent application publication with 
color drawings will be provided by the Office upon request and payment of the 
necessary fee. 

[0017] The present inventive method and its advantages will be better 
understood by referring to the following detailed description and the attached 
drawings in which: 

[0018] Figure 1A and 1B illustrate how Sneirs law of optics can be applied 
to an inclined water bottom to show that a region with lateral velocity variation 
will change the horizontal component of the wave vector of an incident wave; 

[0019] Figure 2 illustrates a velocity model in which the velocity varies with 
both horizontal and vertical position; 

[0020] Figure 3 illustrates an impulse response test at three horizontal 
locations using the velocity model of Fig. 2 and single-velocity pre-stack time 
migration; 

[0021] Figure 4 illustrates the same impulse response test as Fig. 3 but 
with pre-stack time migration corrected by the present inventive method; 

[0022] Figure 5 is a flow chart illustrating the basic steps of the present 
inventive method for correcting pre-stack time migration to treat mild lateral 
variations in subsurface velocity; 

[0023] Figure 6 illustrates the image of a common offset section for a 
point-diffractor model consisting of three point-diffractors at different horizontal 
locations, using pre-stack time migration and a single velocity function 
evaluated at the middle diffractor location; 
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[0024] Figure 7 illustrates the same test of focusing three point-diffractors 
as shown in Fig. 6, except the present inventive method is used to correct 
pre-stack time migration; 

[0025] Figure 8A illustrates a more complex velocity model than Fig. 2; 

[0026] Figure 8B illustrates a layered density model for the same spatial 
region covered by Fig. 8A; 

[0027] Figure 9A illustrates how well pre-stack time migration images 
subsurface reflectors using reflections from the same middle point but with 
different offsets; 

[0028] Figure 9B illustrates the same imaging exercise as in Fig. 9A but 
with the present inventive method used' to correct the pre-stack time 
migration; 

[0029] Figure 10 illustrates the stacked image obtained from single velocity 
pre-stack time migration of the synthetic model data set represented by Fig's. 
8A and 8B; 

[0030] Figure 1 1 illustrates the stacked image of the same synthetic model 
as in Fig. 10 but with the present inventive method used to correct pre-stack 
time migration; and 

[0031] Figure 12 is a flow chart illustrating a process for performing 
controlled amplitude pre-stack time migration in three-dimensional space (the 
"Finn-Winbow" method). 

[0032] The invention will be described in connection with its preferred 
embodiments. However, to the extent that the following detailed description is 
specific to a particular embodiment or a particular use of the invention, this is 
intended to be illustrative only, and is not to be construed as limiting the scope 
of the invention. On the contrary, it is intended to cover all alternatives. 
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modifications, and equivalents which are included within the spirit and scope 
of the invention, as defined by the appended claims. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0033] The present inventive method for correcting mild lateral velocity 
variations in time migration uses several travel time maps obtained using 
velocity functions at locations across a survey to get the phase shift used in 
the migration and the phase shift slope along the lateral direction for a plane 
wave at every vertical time. Both the phase shift and the phase shift slope 
information at each vertical time are used in the migration operator to do the 
migration and correct the lateral velocity variation effect simultaneously in the 
frequency-wavenumber domain. 

[0034] A convenient starting point for explaining the present invention is 
the paper "3-D prestack migration of common-offset data in the frequency- 
wavenumber domain" by Finn and Winbow, SEG International Meeting 
and Exposition, Extended Abstracts (2002). For simplicity, the derivation 
shown below is in two-dimensional space (lateral position x, vertical 
dimension represented by time t). The mathematics for extending to three 
dimensions will be obvious to persons of ordinary skill, and it is intended to 
include three-dimensional applications in this description and in the appended 
claims. 

[0035] Finn's Eq. (4) can be rewritten (in two-dimensional form) as follows: 



where // (kco) is the double Fourier transform of the common offset seismic 
data ju(xj) , CO T (p) is the phase shift for a plane wave with ray parameter 
p = k/co at its stationary point, and y/{x) is the migrated image at horizontal 
position X (and vertical time to, which dependence is implicit in r). The variable 
CD represents the frequency of the wave multiplied by 27r and k is the 
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horizontal (x) component of the wave vector. Finn's amplitude-preserving 
weight function can be and is ignored for purposes of the present invention. 
Finn's integral over o) has been replaced by a finite summation. 

[0036] Next, the migrated image in the ^-domain is obtained by Fourier 
transformation of Eq. (1): 



after switching k and k\ 

[0037] Finn assumed that velocity varies only in the vertical direction which 
means that t(p) depends on the ray parameter p and the vertical time to (and 
not on x). In that event, the integral over x in Eq. (2) above reduces to the 
delta function S{k-k'), which further reduces Eq. (2) to 



Eq. (3) can be referred to as traditional PSTM, as improved by Finn and 
Winbow to accommodate non-zero offsets. 

[0038] The present Invention modifies the approach in Finn to allow mild 
lateral velocity variations to be treated in approximate fashion. Thus, r(p) 
becomes i:(p,x) and Eq. (2) above becomes 




(2) 



(3) 



CO 




(4) 



which can be written 



(5) 



where 
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G{k,k') = p"^(^'^>e'<*'-*'>'dx (6) 
(The integral over A:' was converted to a finite sum.) 

[0039] For general lateral velocity variations, solving Eq. (5) (by computer 
algorithm) is quite expensive due to the double summation. In many real 
situations, however, the velocity has a global, linear increasing or decreasing 
behavior only in a specific lateral direction. In such cases, the phase shift 
o)t(p,x) at vertical time to may be well approximated as a linear function of the 
lateral direction coordinate x, 

T(p, x) = To (p) + 27, (p) X (7 ) 

where coto(p) is the phase shift as a function of ray parameter at a reference 
location, and corj(p) is the phase shift slope along the lateral direction. The 
migration Eq. (5) can then be simplified, beginning with Eq. (6): 

using the definition of the delta function. This reduces Eq. (5) to 

<^W = Z e''''^''''Mik\co) (8) 

where 

k = k'-coT.ip') (9) 

and 

p'^k'/o) (10) 

[0040] By comparing Eq. (8) with Eq. (3), it can be seen that the present 
invention provides a migration operator that corrects for lateral velocity 
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variation by shifting the wave number of the surface seismic wave to a 
different wave number in the image. The amount of shift, cotj (p'J, is the phase 
shift slope which depends on the lateral velocity variation. 

[0041] The derivation above is for any vertical time to. The correction to 
the lateral velocity variation varies with vertical time. In the case of no lateral 
velocity variation, the shift amount of the wave number is zero, and the wave 
number is conserved at all vertical times. The migration operation of Eq. (8) 
becomes Eq. (3). Figures 1A and 1B illustrate the effect of lateral velocity 
variation on the wave vector refracting through a horizontal flat water bottom 
(Fig. 1A) and a dipping water bottom (Fig. 1B). In Fig. 1A, the incident wave 
vector 1 makes an angle with respect to the normal 4 to the horizontal 
water bottom 3. After passing through the water bottom into the subsurface, 
the refracted wave vector 2 now makes a different angle ^ j with respect to 

the perpendicular direction because the wave velocity below the interface is 
different than that above the interface. Snell's Law relates the variables and 
can be written 



where ?^ and ki are the wavelength and wave number (magnitude of the 
wave vector) in water and and k2 are the wavelength and wave number in 
the subsurface below the water bottom. 

[0042] Thus, 



which means that the horizontal component 5 of the wave vector 1 above a 
horizontal water bottom 3 Is equal to the horizontal component 6 of the wave 
vector 2 below the horizontal water bottom. Thus, in the terminology of Eq, 
(9), k = k\ and there is no wave number shift. 



(11) 



sin ^2 ^1 



ki sin^/ = ^2sin<S2 



(12) 
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[0043] Clearly this will not be true for Fig. 1B because Snell's Law applies 
relative to the normal to the interface 8 which is no longer the vertical. In Fig. 
1A, there is no lateral (horizontal) velocity variation; in Fig. 1B, there is a 
lateral velocity variation that occurs where a horizontal line such as 7 
intersects the dipping water bottom 8. 

[0044] In reality, a dipping water bottom typically produces more of a 
lateral (horizontal) variation than merely a step-function change at the water 
bottom. The sediment velocity may often be approximated as a function of 
depth below the water bottom (^-:&^^). Since is a function of horizontal 

position for a dipping water bottom, the velocity throughout the subsurface 
varies continuously with horizontal position at constant ^. Just as in the 
simplified case illustrated in Fig. 1B, the changes in lateral velocity will cause 
changes in the horizontal component of the wave vector, with the general 
result that, for a dipping water bottom, the horizontal wave vector component 
k in the water is not equal to the horizontal wave vector component in the 
sediment. According to the present invention, this wave number shift can be 
approximated by equation (9) for mild lateral velocity variation. 

[0045] In the traditional single velocity function PSTM, the difference 
between k' and k is ignored. In the present inventive method, the migration 
operation of Eq. (8) takes this A: -shift of Eq. (9) into account by calculating the 
phase shift slope as in Eq. (7). By applying this technique, errors in both 
lateral position and vertical time can be reduced greatly. By comparing Eq. 
(8) to Eq. (3), one can see that the increase in computational costs compared 
to the one velocity case is very small. The extra computation is the mapping 
from the wave vector of the surface seismic data to that in the subsurface 
image as described by Eq. (9). 

[0046] As Eq. (7) indicates, the present inventive method works for 
situations where the phase shift changes approximately linearly along a lateral 
direction. In practical applications, a dipping water bottom and/or a 
subsurface velocity with a gradually increasing or decreasing trend along a 
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lateral direction may be well approximated by this linear assumption. One can 
select a few velocity functions across the survey to calculate the phase shift 
tables at those velocity locations. A linear fitting of the phase shift for a given 
ray parameter p and at a given vertical time versus the lateral location x can 
be performed to get the phase shift slope. The migration is performed 
according to Eq's. (8) and (9). 

[0047] A more detailed description of the present inventive method is given 
in the flow chart of Fig. 5. 

[0048] At step 100, the user selects seismic velocity functions v{£) at at 
least two lateral (x) locations. These velocity functions will be different 
functions of ^ . If they were the same function, there would be no change of 
velocity with lateral position (at constant ^ ), and hence there would be no 
need to use the present invention. The number and spacing of lateral 
locations where a v{^) is selected are governed by the need for any linear fit 
to accurately represent the linear component of the data being fit. 

[0049] At step 110, common-offset surface seismic data (gathered from 
the region where the velocity functions in step 100 are evaluated) are 
transformed from the space-time (x,t) domain to the wave number - frequency 
(k,a)) domain. In preferred embodiments of the present inventive method, this 
transformation will be a double Fourier transform. The transformation is made 
because it is well known that a time shift in the x-t domain (i.e., migration of 
the seismic image from the perceived depth to the correct depth) Is equivalent 
to a phase shift in the k-co domain. It is further well known that migration of 
zero offset data by the application of phase shifts in the k-co domain is much 
faster (in computation time) than the equivalent operation of time shifting in 
the X't domain. Finn and Winbow extended this approach to include nonzero 
offset data; however, their method was limited to situations with no velocity 
variation in any lateral direction. The present inventive method relaxes that 
restriction. 
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[0050] At step 120, the user calculates a travel time map for each x- 
location, jc^, for which a velocity function v(^) is selected in step 100. Each 
map determines t{x-xsy^), the travel time from the source location Xs to any 
subterranean point (jc,^). The travel time is calculated by dividing the 
distance traveled (determined from the geometry) by the velocity (assumed to 
be a function only of ^ ). The determination of travel time maps is a well 
known step in performing PSTM. 

[0051] At step 130 of Fig. 5, r-maps corresponding to the travel time 
maps from step 120 are calculated, one such map for each surface (xs) 
location selected in step 100. This is the step discussed previously where a 
correspondence is made between time shifts in the x-t domain and phase 
shifts in the k-co domain. The FInn-Winbow reference mentioned previously 
first disclosed a method for doing this. (Alternatively, see U.S. Patent No. 
6,446,007 to Finn and Winbow.) Their method is restricted to v = v(^), with no 
x dependence, and the present inventive method satisfies that condition since 
a separate r -map based on a one-dimensional velocity function is obtained 
for each surface location Xs, The Finn-Winbow method is described in detail in 
the Appendix hereto. The resulting maps yield r as a function of the ray 
parameter where p^k/w, 

[0052] In each r-map from step 130, r(j)) will not vary with x, but 
between r -maps for adjacent Xs locations, r will vary with jc. (This is why at 
least two Xs locations are required for the present inventive method; a single jc^ 
location is the Finn-Winbow method.) At step 140, r(/?) is found as a linear 
function of jc at each subsurface depth (^) by linear fitting across the r -maps. 
Such linear fits yield the coefficients Tq{p) and T^{p) in Eq. (7) for every 
vertical depth ^ (or vertical time) and every p value. 
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[0053] As previously observed, the migration equation of the present 
inventive method, Eq. (8), is the same as the standard PSTM migration 
equation in the k-co domain (Eq. 3) but with the ^-shift of Eq. (9). With r^ip) 

and Ti (p) now calculated, the difference k - k' can be calculated from Eq. (9) 
and Eq. (8) can then be used to calculate the migrated seismic image in k-o) 
space using well-known PSTM methods. This is step 150 of Fig. 5. 

[0054] In the final step (160), the migrated image is reverse transformed 
(using the transform used in step 110) back to the x-t domain. 

EXAMPLES 

[0055] Three examples are given below to demonstrate the effectiveness 
of the present inventive method for correcting errors in both lateral position 
and vertical time in PSTM in the wave number and frequency domain. Each 
example applies a standard procedure in wide use for testing the accuracy of 
seismic processing algorithms. 

(1) Impulse response test: 

[0056] The first example is an impulse response test for a velocity model 
with a dipping water bottom in two-dimensional space. Figure 2 illustrates this 
velocity model. The water bottom 21 has a dip of about 3 degrees. The 
vertical axis is depth s in meters and the horizontal axis is distance x in 
meters. The water depth ranges from 1500m at the left end to 500m at the 
right. The sediment velocity function is assumed to be given by the following 
linear function of depth ^ : 



where is the water depth in meters. Because the water bottom is slightly 
dipping (approximately 3 degrees), depends on horizontal position x, and 
therefore v is actually a function of both ^ and x. Because the dip is slight, the 



v(5) = 1600 + .4(5 




(13) 
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horizontal velocity dependence is mild, and therefore the situation depicted is 
both common and well suited for application of the present inventive method. 
In the sediment below the water bottom in Fig. 2, the shading depicts (see 
scale 22, in m/s) how the seismic wave velocity in the sediment varies 
continuously with both ^ and x, in accordance with Eq. (13). The seismic 
wave velocity in the water above interface 21 is assumed to be 1500 m/s , and 
is constant throughout the water volume. Impulse responses at three 
horizontal locations, 4,000m, 9,000m, and 14,000m, of two-way travel time 
r = 4 seconds and offset 3,000m are computed and shown in Fig. 3 and Fig. 
4. This means that there are three different source-receiver locations, with 
source and receiver always 3,000 meters apart, and the mid-point between 
them located at 4,000m, 9,000m and 14,000m. The purpose of this particular 
impulse response test is to calculate (and display graphically in Figures 3 and 
4) where subsurface reflectors would have to be located to reflect seismic 
waves originating at the stated source positions back to receivers at the 
stated offset (lateral distance from the source) subject to the constraint that 
the total travel time is always 4 seconds. In Figures 3 and 4, horizontal 
distance in meters is plotted on the horizontal axis, and vertical two-way travel 
time in seconds is plotted on the vertical axis. Both Figures 3 and 4 use and 
depend upon the velocity distribution of Fig. 2. Focusing on the left-most of 
the three concave-up reflectors depicted in either Figure 3 or Figure 4, the 
lowest point lies at about 3.6 seconds on the vertical axis. Since vertical two- 
way travel time is plotted on that axis, a value of 3.6 seconds at x = 4,000 
meters means that the low point of the reflector is located at a depth such that 
an acoustic seismic wave requires 1 .8 seconds to travel vertically downward 
from the water surface with velocity as given by the variable distribution of Fig. 
2. That the "depth" is 1 .8 seconds rather than 2.0 seconds is explained by the 
fact that the source is at x = 2,500 meters and the receiver at x = 5,500 meters 
and thus the detected wave reflected at the low point of the reflector travels 
slightly further than the vertical path for a wave where source and receiver are 
both located at x = 4,000 meters. In this manner, if any point on the reflector 
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curve is connected in a ray path to the source at x = 2,500 meters and to the 
detector at x = 5,500 meters, the total travel time is 4.0 seconds, which is the 
constraint given for the calculational exercise. 

[0057] The other two concave-up curves in Figures 3 and 4 represent 
where reflections would have to occur to generate two-way travel times of 4 
seconds for the other two source-detector locations assumed for the problem: 

X = 7,500m (source) and 10,500m (detector); and 
X = 12,500m (source) and 15,500m (detector) 

[0058] Examining more closely the concave-up shapes in Fig. 3, the black 
and white variable density plot 31 is generated by PSTM using a single 
velocity function (Eq. 13) evaluated at the middle location, x = 9,000m, i.e., 

using Eq. (13) with ^ 1,000/w. The black and white variable density plot is 

used to reflect the finite seismic impulse duration in the model data set. Red 
curve 32 represents the exact result. In Fig. 4, the present inventive method 
has been used to generate the variable density plot. The velocity function of 
Eq. (13) was evaluated at three horizontal locations (4,000m, 9,000m, and 
14,000m) to calculate the phase shift and phase shift slope in Eq. 7. The 
exact result curve is now so completely superimposed upon the variable 
density plot that it is almost impossible to distinguish. Note that the exact 
result (the same as curve 32 in Fig. 3) shows that the reflector bottom falls 
slightly deeper as the horizontal location is moved to a larger x-value, a 
consequence of sediment velocities that increase with x for a given depth (see 
Fig. 2). PSTM is unable to show this effect; the three concave-up variable 
density plots 31 in Fig. 3 are identical in shape and location. There are errors 
in both vertical time and lateral position at locations away from x = 9,000w 
where the velocity function was evaluated. As Fig. 4 shows, the present 
inventive method has compensated these errors caused by lateral velocity 
variations within the migration operation. 
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[0059] For this example, evaluating the velocity function at three locations 
gives excellent results, and there is no need to use more than three locations. 
In general, the number of locations needed to get an acceptable linear fit, 
given the degree of linearity that the data actually exhibit, is problem- 
dependent. 

[0060] The "exact" result as shown by curve 32 in Figures 3 and 4 is 
calculated by well-known ray tracing methods using the full (x,^) dependence 
of the velocity function of Eq. (13). A travel time map is created wherein at 
each point (jc,«) the travel time to (a) the source and (b) the receiver is 
calculated using the velocity function of Eq. (13). Then, all points where the 
two travel times sum to 4 seconds are located and connected with a curve, 
which is curve 32 in Figures 3 and 4. 

2. Focusing of Point Diffraction: 

[0061] The velocity model used in this example is the same as that used in 
the impulse response example, and is shown in Fig. 2. There are three point 
diffractors located at horizontal locations of 4,000m, 9,000m, and 14,000m, 
respectively. A common offset seismic section of offset 3Km was collected. 
This common offset section was migrated by using both the regular PSTM 
and the present inventive method. The images 11, 12 and 13 of the three 
point-dlffractors are shown in Fig. 6 (regular PSTM) and Fig. 7 (present 
invention). The velocity function used in the regular PSTM was taken from 
horizontal location 9,000m. For the present inventive method, three velocity 
functions at horizontal locations 4,000m, 9,000m, and 14,000m were used in 
the calculation of phase shift slope in Eq. 7. A good migration algorithm 
should be able to collapse the diffraction energy into well-focused points. The 
diffractors at horizontal locations 4,000m and 14,000m in Fig. 6 from regular 
PSTM are much less focused compared to the image at 9,000m because of 
the failure to account for the lateral velocity variation in the imaging process. 
All three diffractors are focused very well in Fig. 7 where the A:-shift is applied 
in the PSTM as taught in the present inventive method. 



Express Mail Label No. EU959479572US 



Filed on July 3 1,2003 



Attorney Docket No. PM 99.017 



- 19- 



(3) A General Two-Dimensional Synthetic Model : 

[0062] The velocity and density models used for this example are shown in 
Figures 8A and 8B, respectively. This model has a dipping water bottom 
which can be seen at 81 in both drawings beginning at a depth of 
approximately 90 meters on the left and dipping to approximately 390 meters 
on the right, which is a dip of approximately one degree. As indicated by the 
shading (see scale 82, in m/s), the sediment velocity field is smooth. The 
velocity has both a constant vertical and a constant lateral gradient. The 
velocity contour has an effective dip about 3.6 degrees. A thin salt body 83 
was inserted into the sediment. The density is a layered model. The density 
within a layer is constant. (See density scale 84, in gm/cm^.) The seismic 
reflections come from the density contrast interfaces and the salt boundary. 
A 2-D acoustic wave equation modeling code was used to generate the 
seismic data. Both the single velocity function PSTM (Eq. (3)) and the 
present inventive method (Eq.'s (7), (8) and (9)) have been applied to this 
data set.' The velocity function used in the regular single velocity function 
PSTM was taken from near the center of the model. The migrated common 
middle point gathers from both methods at horizontal location 6,250 meters 
are shown in Fig.'s 9A and 98. The migrated gather from the present 
inventive method shown in Fig. 9B is much flatter than that using the 
traditional single velocity function PSTM shown in Fig. 9A. An example is the 
reflector 92, accurately imaged in Fig. 9B at a vertical two-way time of 
approximately 3.1 seconds. The same reflector is imaged by the curve 91 in 
Fig. 9A. This comparison indicates that the migration operator used in the 
present inventive method has corrected the lateral velocity variation effect, 
and has migrated the reflections of different offsets to the same vertical time. 

[0063] The stacked images of this model from both methods are shown in 
Figures 10 and 1 1 , with Fig. 10 showing the traditional PSTM (as improved by 
Finn and Winbow) method and Fig, 11 the present inventive method. Figures 
10 and 11 reflect the layered density model used in the synthetic data 
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generation, as shown in Fig. 8B. The mis-positioned reflectors can be seen 
clearly in Fig. 10, especially at the left side of the salt body, where the 
reflectors are steep. The phase slope in the present inventive method was 
obtained by using three velocity functions. The imaged reflector positions at 
the left side of the model from the present inventive method in Fig. 1 1 are 
much more accurate (closer to the density model of Fig. 8B) than that of the 
regular single velocity PSTM in Fig. 10. This image improvement is the result 
of including the lateral velocity variations in the migration. Small stairs can be 
detected on the density contrast interfaces in the density model. These small 
stairs give rise to diffraction energy in the data set. To focus the diffraction 
energy correctly, both accurate lateral and vertical positioning are needed. 
Improved imaging of these small stairs in Fig. 11 compared to Fig. 10 can be 
seen by inspecting the reflectors at the right side of the model. 
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APPENDIX 

[0064] According to Finn and Winbow, controlled-amplitude prestack time 
nnigration in the co-k domain is implemented using the following migration 
equation: 

D,(ic)= \d^ke-- /rf.^'-(-)D(^,^-)^^ (A-2) 

where Dm{x) is the migrated image at point 3c, D {oj.k) represents the input 
data after Fourier transformation with respect to time and midpoint so as to 
transform the data to the co-k domain, T{fj,d) is the migration time shift 
function (called r in the description of the present invention), and 
w(77,z)/K3d ifi,^) is an amplitude-preserving weight function. 

[0065] Fig. 12 is a flowchart illustrating the principal steps used in 
performing controlled amplitude prestack time migration in the (o-k domain 
according to Eq. (A-2). At step 46, a set of common-offset seismic data 
traces is obtained. Initially, these data traces are in the space-time {§-t) 
domain. Next, at step 48, a seismic velocity function for the subsurface region 
in question is selected. As noted above, a basic assumption of the Finn- 
Winbow method is that the subsurface seismic velocity is a function of depth 
or time, but is laterally invariant. Techniques for selecting an appropriate 
seismic velocity function for use in the present invention are well known to 
persons skilled in the art and, accordingly, will not be described further herein. 

[0066] At step 50, the common-offset data traces obtained at step 46 
are 3-D Fourier transformed from the f -t domain to the w-k domain. In Eq, 
(A-2), the Fourier transformed data traces are denoted by D(ty,k). 

[0067] At step 52, the seismic velocity function selected at step 48 is 
used to compute a time function for use in migrating the transformed data 
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traces in the o-k domain. This time function is denoted by ^(7,^) in Eq. (A- 
2) and is defined by: 



T{?j.z) = t^-pTj (A-3) 



where 



P = 



(A-4) 



is the midpoint ray parameter computed at fixed depth. As Is well known to 
those skilled in the art, the ray parameter p is the horizontal component of 

the slowness vector. As further described below, the time function T{f],^) is 

converted to the p-^ (dip-depth) domain before it is applied to the 

transformed data traces. 

[0068] At step 54, a weighting function for preserving the seismic 
amplitudes of the input data traces is computed. As noted above, the 
weighting function is denoted by v^{fi,^)/Y^% {fj.^) in Eq. (A-2). 

[0069] is the weight function described by Schleicher et al., in "3-D 

true amplitude finite-offset migration", 58 Geophysics, pp. 1112-1126 (1993). 
According to Schleicher et a!., the weight function depends on <f and x, and 
is given by: 

J| ^\ V cos cos |det(A^^^ + A^^^ } ^^^^^ 

where is the surface dip angle of the ray connecting the source s and the 
reflection point 3c, is the surface dip angle of the ray connecting the 
receiver g and the reflection point jc, is the surface compressional wave 
velocity, and the matrices are given by: 
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(A-6a) 



and, 



d% 



d4^dp dtdp 



(A-6b) 



where p and a are coordinates on a plane tangential to the reflector at the 
reflection point R. Given a knowledge (for example, from a time map) of the 
rays connecting the source and receiver to the reflection point, all of the 
above derivatives can be calculated. The source and receiver locations are 
s=i + h and g = ^-h, where h Is the half-offset as a two-dimensional 

vector on the earth's surface. Since ^ = 3c - 1 , the weight function \^,x) may 
be expressed as a function of ^ and 3c or as a function of fj and s. 
Alternatively, the weight function w{fj,^) can be computed using dynamic ray 
tracing, as described In Hanitzsch, "Comparison of weights in prestacl< 
amplitude preserving Kirchhoff depth migration," 62 Geophysics, pp. 1812- 
1816 (1997). 
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[0070] The factor K^j^ifj,^) is a normalization factor used to adjust for 
differences between working in 3c space and working in k space. K^j^ifj,^ is 
defined by the following equation: 



^3i>(^>^) = det 



dTj.dijj 



(A-7) 



where the indices / and J take on the values x and y corresponding to the two 
components of the variable ij . 

[0071] As with the time function T{fj,^), the amplitude-preserving 

weighting function v/{rj,z)/K^^^{f},s) is converted to the p-^ (dip-depth) 

domain before it is applied to the transformed data traces. This conversion is 
described in detail below. 

[0072] At step 56, the transformed data traces D[co,k) are phase 
shifted using the time function calculated at step 52. 

[0073] Next, at step 58, the phase-shifted data traces are multiplied by 

the amplitude-preserving weighting function calculated at step 54. 

[0074] Finally, at steps 60 and 61 , the phase-shifted and weighted data 
traces are transformed from the co-k domain back to the x-^ domain. This 
is accomplished in two stages. First, at step 60, an image is formed in the 
wavenumber {k) domain by summing with respect to frequency co. Then, at 

step 61, the k domain image is 2-D Fourier transformed back to the x 
domain. The result of this process is x-^ domain migrated data traces in 
which the migration process has preserved seismic amplitudes. 

[0075] The functions T and w/kI^J are together referred to as the 
"migration operator," M = {t,w/k1^j^). A fundamental difficulty in the 
implementation of Eq. (A-2) arises because the functions T and w/kI^^ are 
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initially calculated in fj-s space while the transformed data traces are \n co-k 
space. To overcome this difficulty, rand w/k\^^ are first converted to the p - 
z domain, as further described below. The conversion of r and w/k\'^ \o 
5 space permits an efficient implementation of Eq. (A-2) since points in p-^ 
space can be related to points \t\ ay-k space using the following relationship: 

P = - (A-8) 
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